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Abstract 

We present an exact solution for an asymmetric exclusion process 
on a ring with three classes of particles and vacancies. Using a matrix 
product Ansatz, we find explicit expressions for the weights of the 
configurations in the stationary state. The solution involves tensor 
products of quadratic algebras. 

1 Introduction 

The one-dimensional asymmetric simple exclusion process (ASEP) has been 
extensively studied in mathematical and physical literature as one of the 
simplest models for non-equilibrium statistical mechanics |l], ||, [| ||| . The 
ASEP is a model of particles diffusing on a lattice driven by an external 
field and with hard-core exclusion. A variety of different phenomena can be 
described by the exclusion process, for instance superionic conductors ||, 
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traffic flows Q and interface growth Exact results have been obtained 
for the one-dimensional exclusion process with the help of various methods 
such as the Bethe AnsatzjH, ^, pL0f| , and, more recently, a matrix product 
Ansatz (see for example fllTf ), 

The matrix Ansatz has led to new exact results concerning the station- 
ary state of various models. Originally developed for the study of directed 
animals on a lattice [12^1 , this method has been successfully applied to the 
exclusion process with open boundaries |13|]. It has then been extended to 
study systems with second class particles and shocks [|l4|] , time discrete dy- 



namics 1151 16, 17, 18, 19, p(J|, and to calculate diffusion constants plf. The 



algebras involved have led to interesting representation problems j22|. 

Models with more than two classes of particles have hardly been inves- 
tigated (| 23], Here, we study an exclusion model with vacancies and 
three classes of particles on a ring. Up to now, it was not known whether the 
matrix Ansatz could be used to construct an exact solution of this model. In 
this paper, we shall give an exact expression for the stationary state of our 
model by using a suitable matrix Ansatz that involves tensor products of 
quadratic algebras. Apart from the question as to how far the matrix Ansatz 
can be used, the model is interesting in itself, because it is suitable for a 
detailed study of shocks p5f| . In section 2, we define the model and explain 
what the matrix Ansatz is. In section 3 we recall the matrix solution for the 
asymmetric exclusion process with second class particles. In section 4, we 
give an explicit representation of the algebra that describes the stationary 
state of our model and present a proof of our solution. We also give an 
explicit solution for the stationary state without using any representations. 
We then discuss an algorithm that allows to obtain exact properties of the 
stationary state for large systems by using a computer. The concluding sec- 
tion discusses our results and some generalizations. The appendix contains 
details of the proof and certain algebraic properties and recursion relations. 



2 Definition of the model 

We consider a periodic one dimensional lattice of L sites. Each site of the 
lattice is either empty or occupied by one single particle that can be of type 
1,2 or 3. For reasons that will become apparent later on, we say that empty 
sites are occupied by holes (or vacancies) and we shall call holes as particles 
of the fourth type. We denote the number of particles of type k in the system 
by rife, where k = 1,2,3,4. The state of a site i is specified by a discrete 
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variable n that takes four possible values: 

Ti = 1, 2, 3 or 4 if site i is occupied by a particle of type 1, 2, 3 or 4 . (1) 

The dynamics of the system is given by certain transition rules. During 
an infinitesimal time step dt, the following processes take place on a bond 
(i, i + 1) with probability dt: 



12 -> 21 with rate 1 

13 — > 31 with rate 1 

14 — ► 41 with rate 1 



23 32 with rate 1 

24 42 with rate 1 



34 43 with rate 1 . (2) 

All other transitions are forbidden. Obviously, the dynamics conserves the 
number of particles and one has J2k=i n k = L. It is also clear from these 
rules that particles of type n can 'overtake' particles of type type m only 
if n < m. The transition rules therefore induce a hierarchy among the 
particles. 

In the literature, particles of type k = 1,2,3 are named first, second 
and third class particles. The model defined in (g) is called 'the asymmetric 
exclusion process with three classes of particles and holes'. 

Note that a first class particle behaves always in the same way in regards 
to all the other particles, whereas a third class particle for example behaves 
like a first class particle in respect to the holes, but as a hole in respect to 
the second and first class particles. 

The rules given in @ are translationally invariant. Using this property, 
we decide that a particle of the third class occupies the site number L and 
we enumerate the different configurations of the system. The total number 
of configurations N to t is given by 

_ (L-l)! 
tot " m!n 2 !(n 3 - l)!n 4 ! ' [ ) 

The dynamics of the system can be encoded in a Markov matrix M of 
size iVtot x ^tot- T ne coefficient M{C,C) of this matrix represents the rate 
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of transition from a configuration C to a different configuration C'; M(C,C) 
is the rate of exit from a given configuration C. Due to the local structure 
of the rules (0) , M can be written as a sum of local operators that represent 
the transitions that take place at a bond (i, i + 1) 

L 

M = ^m M+ i . (4) 

i=l 

An explicit expression of the matrices m^j+i is given in the appendix. 

In the long time limit, the process reaches a stationary state in which 
each configuration C of the system has a stationary probability p(C). The 
stationary state exists and is unique. This follows from the fact that the 
rules @ define an irreducible Markov process, i.e. any given configuration 
can evolve to any other one. The properties of the stationary state can be 
determined once the probabilities p(C) are known for all C. In equilibrium 
statistical mechanics these numbers are given by the Boltzmann factor, but 
in our model there is a priori no method to calculate these quantities: one 
has to solve the stationary master equation 

J2M(C,C')p(C') = . (5) 

c 

This is a system of N^ Q ^ coupled linear equations whose complexity grows 
exponentially with the size L of the system. 

The matrix Ansatz consists in solving the system @ by writing the 
probabilities p(C) as a trace of a product of four non-commuting operators, 
Ai, A2, A3 and A4, where each represents one type of particle: 

p(C) = ^Tr(A T1 ...A TL ) (6) 

where A Ti is equal to A^ if site i is occupied by a particle of type k (k = 
1,2,3,4) in configuration C. The constant Z is an overall normalization 
factor, that depends on L and on the n^'s; it ensures that J2p(C) = 1. 

In section 4, we shall present explicit operators A^ {k = 1,2,3,4) and 
prove that the weights p(C) constructed from @ are solutions of the master 
equation (||). Certain properties of the solution of the model with just first 
and second class particles (i.e in our case 713 = ) are useful for constructing 
the operators A^. Therefore, we review this solution in the following section. 
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3 Matrix solution of the ASEP with first and sec- 
ond class particles 

Suppose first that there are only first class particles and holes on the ring. 
In this case, the steady state is such that all configurations C have the 
same weight |26fl . Assuming that the last site is always occupied by a first 
class particle, we obtain that the total number of configurations will be 
(L-nf)7(nl-i)! anc ^ eacn configuration has a stationary probability equal to: 



P(C) 



(m-l)!(L-m)! 
(L-l)! 



(7) 



In this simple case, a matrix Ansatz is not needed (one can choose the 
matrices representing particles and holes to be both equal to the scalar 1). 

We consider now the model defined in (|2|) without particles of type 3 
(i.e. 713 = 0). There are n\ particles of the first class and n-i particles of 
the second class. For this model, the stationary probability is non-uniform 
and was obtained in [14| from a matrix product Ansatz. Following [14], we 
denote by D, E and A the operators that represent particles of type 1, 2 and 
holes, respectively. The numbers p(C) obtained from expression (|J) are the 
stationary probabilities of the exclusion process with first and second class 
particles, if the three operators E, D and A generate the quadratic algebra 
defined by the relations 



BE 
DA 
AE 



D + E 

A 

E 



(8) 



It is convenient to work with an explicit representation of the algebra 
(H). A particularly useful choice is: 
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A=|l)(l| 



( 1 







V • 



(9) 



The operators D and E are represented by matrices that act on an 
infinite dimensional space with denumerable basis (|1), |2), ...|n), ...). The 
operator A is a projector of rank 1 on the first element of the basis and has 
a finite trace. This ensures that any expression of the type (|6|) is finite. 

Using the algebraic rules (|8|) or the explicit representation (|9|) , all station- 
ary probabilities are determined. In order to calculate physical quantities 
such as density profiles, or average local currents, we must know the con- 
stant Z, which plays a role analogous to that of the partition function. The 
expression for Z is simple if there is only one second class particle in the 
system, always located on the last site. In that case, one obtains 



Z 



When the density of second class particles is finite, asymptotic formulae for 
Z are derived for systems of large size, using a grand canonical formalism 
Ll. 




4 Matrix solution of the ASEP with first, second 
and third class particles 

4.1 Explicit forms of the matrices 

We shall show that the stationary weights, solutions of the master equation 
(||), can be calculated via the Ansatz @, from the following four operators 
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All matrices are infinite dimensional and their coefficients are themselves the 
infinite dimensional operators D and E of (||) which satisfy DE = D + E 
and do not commute with each other (i.e. scalar representations of D and 
E are excluded). Another way to look at the operators given in ( |i~0|) is 
to consider them as matrices operating on an infinite dimensional space, 
with non-commutative elements. The operators A 2 and A3 have only two 
non-zero elements and the following relation holds: 



AoA 



2^3 



( A 









v 



fill 



/ 



here A is the rank one projector of @. Before we prove that the stationary 
probabilities given in terms of the A^s solve the master equation (|5|), we 
have to ensure that they are finite. This is not obvious because none of the 
operators given in (fill) has a finite trace. 



4.2 Proof of the finiteness of the Ansatz 

We rewrite the expression (|6|) for the stationary weights as follows 

p(C) = ±Tr(A Tl ...A TL ) = ±Tr(YA 2 XA 3 T) (12) 

where Y and T are products of the A^s (k = 1, 2, 3, 4) and X is a product 
of p (p < L — 2) matrices A± and A^ only. Such a factorization is possible: 
the term A2XA3 starts from the furthermost (proceeding from left to right) 
factor A2 in (A T1 ...A TL ) and ends when an A3 appears for the first time after 
this A 2 - Such a factor A3 always exists since A TL = A3. 

We now prove that the matrix A2XA3, where X is a product of p factors 
A\ and A4, can have non-zero elements only in its first (p + 2) lines or 
columns. We shall say, in such a case, that A2XA3 is 'of finite size (p + 2)'. 
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The operators Aq and A± have two invariant subspaces, the subspace 
generated by the odd vectors of the basis (|1), |3), ...\2n + 1),..-) and the 
subspace generated by (|2), |4), ...|2n), ...). The action of A\ (and that of 
A4) on both invariant subspaces is the same. Therefore, the product X will 
be represented by the following matrix: 



X 
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• / 



(13) 



The symbol * denotes unspecified matrix elements. We emphasize that the 
coefficients (1,1) and (2,2) of X are identical. This coefficient is a matrix 
X which is a linear combination of various products, each product having p 
factors, and each factor being either a D or an E. 

Using the expressions of A2 and ^3 and the matrix (O) for X, we find 



A 2 XA 3 



( D X E - E X D 









V 



(14) 



the product A2XA3 has only one non-zero coefficient D\E — E\D where 
X is a linear combination of products of p factors D and E. Therefore, we 
need to show that if M is any product of p factors D and E, the matrix 
DME — EMD is finite of size (p + 2) at the most. This is achieved by 
induction on p and by using the explicit representations of D and E given 
in (|). 

For p = 0, M = 1 we obtain DE — ED = A which is a matrix of size 1. 
Now suppose that that our assertion is true for (p — 1). Then, let the matrix 
M be a product of p factors D and E; if M = DM\ (the case M = EM\ is 
similar), we have 

DME - EMD = D{DM l E - EM X D) + (DE - ED)M\D . 

By the induction hypothesis, DM\E — EM\D is finite of size (p + 1) and 
multiplying it by D will increase its size by at most 1 (one verifies this by 



S 



using the explicit representation given in @). The operator (DE—ED)M\D 
is equal to AM\D and is of size less than or equal to (p + 2). 

We have shown that the factor A2XA3 in (12) is of finite size. Multi- 
plying it on the left or on the right by any of the operators A%,A2, A3 and 
A± given by ([ji]) does not alter this property since the A^s are made of D 
and E's and have only a finite number of non-zero coefficients in each line 
and each column. This proves that the matrix (A T1 ...A TL ) has a finite trace 
and that the stationary probabilities given by (pi) are well defined. 



4.3 Proof of the Ansatz 



We shall use the technique developed in |28| (see for example [29] for details). 
We represent the collection of the (unnormalized) stationary weights p(C) 
as a state vector 

\p) = Tr(^ L ), (15) 

where <S> denotes the tensor product and A is a column vector, having ma- 
trices as entries: 

/ M \ 

A 2 
A3 



A 



(16) 



This allows us to interpret the Markov equation (||) as a stationary 
Schrodinger equation with the non-hermitian 'Hamiltonian' M 



M\p) = J2m,i+l\p) = 



(17) 



i=l 



The matrix m^j+i is local and acts only on the ith and the (i + l)th element 
of the tensor product in (|l5|). We show that m^j+if^l®^] is a divergence-like 
term, i.e. there exists a vector A defined analogously to A 



A 



( i i \ 

A 2 

A3 

V M J 



such that: 



mii+i [A®A] = A<8)A — A<8)A . 



(18) 



(19) 



Summation over i leads to a global cancellation, proving thereby that the 
Markov equation (|l7]) is satisfied. The proof rests upon finding four matrices 
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Aq, A\, A2 and ^3 that satisfy equation (|l9|). In the appendix, we write the 
16 quadratic equations that couple the ^4 K 's and the A^s (see equations 
A.1-A.3). An explicit representation of the A^s that solves these equations 
is given below (here 1 denotes the identity matrix) 
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(20) 



/ 



Remark: There is one subtlety involved here. One should not verify the 
cancellation mechanism only formally but also make sure that all the traces 
of all the products in (|l9|) exist. In all cases but one, this follows from section 
4.2. However, the case when the last factor of the trace is ^2^3 needs extra 
care because the dynamics permutes these two factors. One needs to prove 
the following relation, where Y denotes any product of the matrices A^, 
k = 1,2,3,4: 



Tr{y(^ 3 i 2 - 4^2)} = Tr{A 2 YA 3 - A 2 YA 3 )} 



(21) 



However, this identity can be proved via a reasoning similar to that of section 
4.2. 
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4.4 Representation— free solution 



Consider the following choice for the Ai and Ai operators: 



A 2 
A 3 
At 



1 ® D + (D - l) 2 ® £ , 
A ® D + - 1)) ® , 

1 ® D + (£ - l) 2 ® D 



(22) 
(23) 
(24) 
(25) 



and 



Ax 
A 2 
A 3 
M 



1 <g> (1 + D/2) + (D - l) 2 8) (D/2 - 1) 
A® (1 - D/2) + (iA(D - 1)) ® (£7/2 - 
A ® (£7/2 - 1) + (i(l - £)A) ® (D/2 - 
-{1 ® (1 + E/2) + (E — l) 2 ® (D/2 - 



-1) 
1) 
1)} 



(26) 
(27) 
(28) 
(29) 



where again DE = E + D, A = DE — ED (projector) and where i = 
It is then possible to show that these operators solve (A.1-A.3) and that, 
assuming a third class particles at site L, all the traces in @ and (|l9|) 
are real and finite. The calculation is very similar to the one which was 
presented in the preceding subsection and will therefore be omitted here. 
We believe that this representation-free solution will help to generalize our 
solution to the case of a model with N types of particles. 

4.5 Algebraic properties and recursion relations 

From our solution it is straightforward to derive certain algebraic properties 
of the Ai operators and therefore to find recursion relations of the relative 
weights of the configurations in the system size. Some of those relations 
are listed in appendix B. In fact, we found the operators Ai by solving the 
model for small system sizes on the computer, guessing recursion relations 
and constructing suitable operators which fulfilled these relations. We want 
to remark that it seems to us extremely unlikely that a solution could have 
been constructed just by inspection of equations (A.1-A.3) of Appendix A. 
However, these equations turned out to be very useful for proving that the 
weights given in terms of the A^s are indeed the stationary weights. 

4.6 Finite size cut-off 

We have shown that the matrix Ansatz using the operators given in (|l^) is 
well defined and satisfies the master equation. How can the representation 
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be used for actual computations for systems of size LI The proof of the 
finiteness of the trace (section 4.2) provides a method to compute numer- 
ically the weight of any configuration of size L without involving infinite 
matrices. We showed that all the matrices used to evaluate the trace are at 
most of size L. Therefore, the operators A^s, D and E can be restricted to 
a finite size A, with A large enough to ensure that the Lx L matrices needed 
to calculate the weights are the same as those obtained by multiplying in- 
finite dimensional matrices. Such a cut-off procedure is possible due to the 
bidiagonal structure of the ^4 K 's and of D and E. For example if we limit 
D and £ to a finite size N and consider a product of p such matrices, the 
(N — p — 1) first rows and columns of the resultant matrix will be the same 
as those obtained by taking the product of the initial infinite dimensional 
matrices. To compute the weights of systems of size less than L, we must 
take A > 2L. The computation time using the matrix Ansatz grows alge- 
braically with the system size whereas the increase is exponential for solving 
the master equation. Thus our solution allows an exact numerical study of 
such systems for large sizes ]25[| . 

To determine the currents and the density profiles, one needs the normal- 
ization factor Z. Using exact results for systems of sizes up to 8, we guessed 
the following formula for Z for the case when there is only one particle of 
the third type (n^ = 1) and 711,77,2,774 ^ 




5 Conclusion 

We have studied a generalization of the asymmetric exclusion process to a 
system with three classes of particles and holes. This model can be mapped 
to an integrable two dimensional vertex model of equilibrium statistical 
physics |3(| , but the Bethe Ansatz does not allow a simple determination of 
the ground state of this vertex model. However, the stationary weights can 
be calculated via a matrix product Ansatz. Although analytical formulae 
may be difficult to derive (the computation of the constant Z will require 
calculations similar to those of the diffusion constant of an open system |^7|| ) 
the matrix Ansatz enables a much faster exact numerical computation of the 
stationary state of finite size systems. 

Our main interest is theoretical. We wanted to know what kind of al- 
gebras (if any) appear in multi-species processes. Some authors 23, [31], 



32 1 have used generalized quadratic algebras to study systems with many 
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species. Associativity [23| and finite trace condition |3^] for these algebras 
impose severe restrictions on the rates of exchange between different types 
of particles. The simple rates we choose (|2|) do not satisfy these limitations. 
Hence the algebra we have found is not quadratic but rather a tensor prod- 
uct of quadratic algebras. As emphasized in section 4.5 and in Appendix B, 
the identities that are satisfied by the matrices A R can be cubic, quartic or 
of any higher order. We believe that the tensor structure we have obtained 
is fairly general. If for some special choices of the transition rates in (0), 



the matrices D and E can be taken to be scalars [13], our matrices ( |i0|) will 
generate a quadratic algebra. 

There is a recursive structure when one adds new types of particles. The 
exclusion process with only one class of particles is solved by taking the 
matrices representing holes and particles to be both equal to 1. For two 
classes, the matrices D, E and A are infinite dimensional matrices with l's 
as coefficients. For the three classes problem, the matrices given in (|i"o| ) are 
infinite dimensional with D and E as coefficients. 

It is therefore natural to define a generalization of the exclusion process 
for N types of particles, with a priority rule such that a particle of type 
n can overtake a particle of type m if and only if n < m. This model is 
still integrable, and some exact results can be obtained via a Bethe Ansatz. 
Besides, from numerical solution of small systems one finds many relations 



between the rates and the matrices representing each type of particle [33]. 
We hope that our solution will help to find a solution for this generalized 
problem. 

We have studied only the totally asymmetric exclusion process, it may be 
interesting to try to solve the partially asymmetric exclusion process where 
all the rules in (0), such as '12 — > 21 with rate 1', are modified as follows 

12 — > 21 with rate p 

21 -> 12 with rate q (30) 

with p + q = 1. We believe that a suitable tensor product structure should 
allow to compute the ground state of this model. Since such a model could 
presumably display spontaneous symmetry breaking [34|, this would be of 
special interest. 
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A Explicit form of the local Markov matrix 

The local Markov matrix m^j+i that describes the updating of a bond + 
1) is given in the canonical basis (11), (12), (13), (14), (21), (22), .., (44) by a 
16 x 16 matrix: 



mi,i+l = 
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The local divergence condition (|T^) translates into the following 16 cou- 
pled quadratic equations: 

AiAj = AiAj-AiAj for i<j, (A.l) 

AiAj = AjAi - AjAi for i < j , (A. 2) 

= AiAi - AiAi for all i , (A.3) 

where 1 < i, j < 4. 
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B Algebraic properties and recursion relations 



The matrix algebra method is a way to encode recursion relations between 
stationary probabilities of systems of size L and of systems of size (L — 1). 
For some simple models, the matrices can be constructed using 'empirical' 
recursion relations observed on exact solutions for small systems. In our 
case, a complete set of such relations between size L and size (L — 1) is 
difficult to obtain. However the matrices A X ,A 2 ,A 3 and A 4 given in ( |l0| ) 
satisfy a number of algebraic identities that allow to deduce some recursions 
between system of different sizes. We now describe some of these relations 
that generalize the simple quadratic algebra (||). 

1. The matrices A\ and A4 satisfy the algebra that was found in [27| and 
was used to compute the diffusion constant for an open system: 

A 1 A p f 1 (A A A l - AiA 4 )A q 1 ~ 1 A 4 = 
A P 4 A\A 4 - A X A P ± X A\A± - A^A^ 1 A 4 + A X A\A\ (B.l) 

where p and q are strictly positive integers. 

2. Some relations reduce the system size and are reminiscent of the DE = 
D + E identity in (||). However, the following relations are cubic and not 
quadratic. 

^2^2-44 = -42-42 + A 2 A 4 

A X A 3 A 3 = A X A 3 + A 3 A 3 (B.2) 

3. Other relations are similar to the second and the third equality in (^|): 

A 2 A 2 A 3 = A 2 A 3 
A 2 A 3 A 3 = A 2 A 3 
A 2 A 3 A 2 A 3 = A 2 A 3 (B.3) 

This last equality shows that the operator (^2-43) is a projector as we noted 
in (0). 

4. Some rules transform some particles into others without reducing the size 
of the system: 

= -4 2 -4 2 
{A x Ai -A i A 1 )A 2 = A 2 A 4 A 2 

A 3 Ai = A 3 A 3 
A 3 (A X A A - A A A X ) = A 3 A X A 3 

A 3 A 2 A 4 = ^3^2-43 + ^3^3^2 

A X A 3 A 2 = A 2 A 3 A 2 + A 3 A 2 A 2 (B.4) 
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